Draft version September 4, 2008 

Preprint typeset using I4TgX style emulateapj v. 1 0/09/06 



00 

O 
O 
(N 

D 



THE DYNAMICAL STRUCTURE OF DARK MATTER HALOS WITH UNIVERSAL PROPERTIES 

Emmanuel Van Hese, Maarten Baes and Herwig Dejonghe 

Sterrenkundig Observatorium. Universiteit Gent, Krijgslaan 281 S9, B-9000 Gent, Belgium 
Draft version September 4, 2008 

ABSTRACT 

A^-body simulations have unveiled several apparently universal properties of dark matter halos, including 
a cusped density profile, a power-law pseudo phase-space density p/cr^, and a linear /3-7 relation between 
the density slope and the velocity anisotropy. We present a family of self-consistent phase-space distribution 
functions F{E,L), based on the Dehnen-McLaughlin Jeans models, that incorporate these universal properties 
very accurately. These distribution functions, derived using a quadratic programming technique, are analytical, 
positive and smooth over the entire phase space and are able to generate four-parameter velocity anisotropy 
profiles (3(r) with arbitrary asymptotic values /3o and /3oo- We discuss the orbital structure of six radially 
anisotropic systems in detail and argue that, apart from its use for generating initial conditions for V-body 
studies, our dynamical modeling provides a valuable complementary approach to understand the processes 
involved in the formation of dark matter halos. 

Subject headings: cosmology: dark matter - galaxies: clusters: general - galaxies: kinematics and dynamics - 
methods: analytical 
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1. INTRODUCTION 

In theoretical astrophysics the steady increase in computa- 
tional power has sparked a proportional interest and progress 
in the study of large-scale structure formation. In particular, 
as V-body simulations of cold dark matter halos have become 
more detailed, several "universal" properties have emerged. 
We highlight three important characteristics. 

Firstly, numerous cosmological stud ies (e.g. 
Dubinski & Carlber2 1991; Crone et al. 1994; Navarro et al] 
]_996; Fukushige & Makino 1997; Navarro et al. 1997; 
Carlberg et al.. ,1997; .Moore et ah ,1998, .1999: .Jing & Sutg 
20001) revealed similar density profiles over several orders 



of magnitude in halo mass, with a central cusp and an a 
p(r) cx falloff at large radii. These generalized NEW 
models can be described by a general 3 -parame ter family, 
referred to as the Zhao models (or a/37-models) (iHernquistI 
119901; IZhaolll996h . They are defined by the density 



2(700-70)/!; , 



p(r)= 

(rA.)^°(l+(rA. 

or in terms of the logarithmic slope. 



, >)\(7oo-7o)/i; ' 



7(r) = 



dlnp^ ^ 7o + 7< 

in = 

dlnr 



l + ('-A.s)' 



(1) 



(2) 



Although in recent years alternative profiles base d on the Ser- 
sic law have produced e qually good results ( Navarro et alJ 
l2004t iMerritt et al.ll2005l) . the generalized NEW models re- 
main very popular and successful to represent dark matter ha- 
los. 

A second relation was found by Tavl or & Navarrol (l2001h . 
These authors identified that the quantity Q{r) = p/cr'ir), 
which has become known as the pseudo phase-space density, 
behaves as a power law over 2-3 orders of magnitude in radius 
inside the virial radius. 



Q(r) oc r 



(3) 
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Other studies (e.g. iRasia et al.l 120041; lAscasibar et al.l |2004 
have confirmed the scale-free nature of Qif), and their results 
indicate that its slope lies in the range a = 1.90 ±0.05. This 
property is remarkable since the density p{r) nor the velocity 
dispersion air) separately show a power-law behavior 

Finally, the velocity anisotropy profiles /3(r) of dark matter 
systems also evolve toward a similar shape, steepening grad- 
ually from isot ropic in the center to radi ally anisotropic in the 
outer regions. [Hansen & Moord (l2006h suggest a nearly lin- 
ear relation between the logarithmic density slope 7(r) and the 
velocity anisotropy profile, based on various types of equili- 
brated simulations. They proposed the /3- 7 relation 



/?(7)~ 1-1.15(1+7/6). 



(4) 



Several theoretical studies have been made to investigate 
whether solutions of the Jeans equation exist that encompass 
the observed properties of dark matter halos. In particular, 
iDehnen & McLaughlin (2005) investigated the anisotropic 
Jeans equation constrained by a slightly different form of the 
pseudo phase-space density, namely Qrir) = pjo]. with (T,(r) 
the radial velocity dispersion. They found a special solution, 
namely an analytical self-consistent potential-density pair of 
the form ([T]| with an anisotropy profile 

(5, 

1 + ('"/ '"a) 

that also has an exactly linear /3 - 7 relation (in other words, 
= Ts and 2(5 = r\). As these Jeans models satisfy the three 
universal relations mentioned above, we can consider them as 
representative models for realistic dark matter halos. 

The goal of this paper is to take the analytical study of dark 
matter systems a step further, i.e. to look for full dynamical 
models that encompass the universal properties found in V- 
body simulations. In concreto we will look for phase-space 
distribution functions (DFs) v) that self-consistently gen- 
erate the required density, potential, and anisotropy profiles 
encountered in dark matter halos. Such dynamical models 
would provide a very useful complementary approach to gain 
insight into the structure of dark matter halos. 

We shall focus on spherical models, for which a dynami- 
cal description simplifies substantially as in this case the DF 
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can be expressed as a function F{E,L) of the binding energy 
and the angular momentum. But even if we limit ourselves 
to the spherical case, it is not straightforward to obtain such 
dynamical models. A simple approach is to solve the Jeans 
equation and approximate the velocity dispersion pro files by 
a multivariate Gaussian ( Hernquist 1993) . Kazantzidis et alJ 
(|2004) demonstrated however that these systems are far from 
equilibrium, thus losing their initial dynamical structure as 
they evolve to non-Gaussian velocity distributions in subse- 
quent A^-body simulations. We therefore need the tools to de- 
rive full self-consistent equilibrium models, governed by DFs 
with a sufficiently general velocity anisotropy profile. 

Simple analytical models are found for only a few spe- 
cial cases, such as the Plummer, He r nquist, iso c hrone , 
and 7-models ( e.g. Dejonghe 1987; Hernquist 19901 
Baes & Deion ghd 12002; An & Evans 2006; Buyle et all 



2007ah . none of which are able to describe reaUstic dark mat- 



ter halos. For gener al potential-density pairs. lOsipkovl ( Il979h 
and iMerrittj (Il985h provided an algorithm to construct dy- 
namical models with a specific velocity anisotropy profile of 
the form (|5]l, with (3q = 0, /3oo = 1 and 6=1, leaving r.^ as a 
free parameter. Whil e the Osipkov-Merritt method has been 
widelv adopted (e.g. I^r owll2000t iL okas & Mamonl l2001t 
iKazantzidis et al.l|2004l) . comparison with simulations shows 
that such anisotropy profiles are too steep to describe dark 
matter (Mamon & Lokas 2005). Halos are not completely ra- 
dial at infinity (i.e. /3oo < 1) and the linear /3-7 relation can 
only be realized with a lower transition rate 6 <\. On a dy- 
namical note, the associated DFs are of the form F{Q) with 
Q = E-l? I2r\, creating an unphysical cut-off boundary for 
orbits with 2 < 0. Hence, the Osipkov-Merritt framework is 
too limited to generate realistic dark matter systems, and a 
more extens ive method is needed. 

Recently, IWojtak et al.l (l2008l) presented an interesting ap- 
proach to generate dynamical models for potential-density 
pairs with a more general anisotropy profile. They pro- 
posed to express the DF as a separable function of the form 
/£(£')/l(L) where /L(i) is a double power-law function with 
three parameters /3o, /3oo, and Lq- Once their values have 
been determined, the function fsiE) is derived from the ob- 
served density profile by a numerical inversion. This tech- 
nique yields a three-parameter anisotropy profile that resem- 
bles eq. (|5]i, where Lq has a similar role as r^, and a fixed 
transition rate 0.5 < (5 < 1. In this manner, the authors were 
able to construct models with an NEW density with velocity 
dispersion profiles that agreed with their dark matter simula- 
tions. 

In this paper, we present a technique that enables us to 
obtain dynamical models with exactly the four-parameter 
anis otropy profile (|5]l. We demonstrated in a previous pa- 
per (iBaes & Van Hese 2007, hereafter Paper 1) how this can 
be achieved, by postulating a separable parameterized form 
of the so-called augmented density, which provides an equiv- 
alent description of a dynamical system. In certain special 
cases the transformation from the augmented density to the 
DF is analytically tractable. We were thus able to derive 
a family of DFs for the generalized Plummer models (also 
called a-models or Veltmann models, .Veltmann (,1979) ) with 
a linear /3-7 relation. 

Since the generalized Plummer potential-density pairs form 
a special class of the Zhao models ([T]), we will now demon- 
strate that this result is a first step toward more representative 
dark matter profiles. In particular, we seek a family of DFs 
for the Dehnen-McLaughlin systems, because of their unique 



property to satisfy a universal density, a power law Q, {r), and 
a linea r 0-^ relation. W e use a quadratic programming tech- 
nique (iDeionghdl 1989 1) to build DFs as a linear combination 
of base functions of the form derived in Paper I. In this man- 
ner, we demonstrate that it is indeed possible to generate full 
dynamical models for the Dehnen-McLaughlin Jeans models, 
thus providing DFs that encompass the observed properties of 
dark matter halos. 

Our paper is organized as follows. In Section |2] we de- 
scribe the notion of a dynamical model and we briefly summa- 
rize the main aspects of the Dehnen-McLaughlin halos. Next 
we outline our modeling technique in Section |3] we explain 
the quadratic programming algorithm and we recapitulate the 
functions that we derived in Paper 1. With these components, 
we can build a library to apply the QP-method to the Dehnen- 
McLaughlin halos. In Section |4] we present our results for a 
set of systems with different velocity anisotropy profiles, and 
we discuss the moments, phase-space DFs, and energy and 
angular momentum distributions for these models. Finally, 
we formulate our conclusions in Section|5] 

2. PRELIMINARIES 

2.1. Dynamical models 

The dynamical structure of a gravitational equilibrium sys- 
tem is completely determined by the DF F{r,v), which 
describes the probability distribution of particles in six- 
dimensional phase space. A dynamical model is only physical 
if its DF is non-negative everywhere. In the case of spherical 
symmetry, this DF can be written as a function F{E,L) of two 
isolating integrals, namely the binding energy and the angular 
momentum 



1,1, 

L=rvT, 



with 



Vj : 



(6) 
(7) 

(8) 



the transverse velocity, and ipir) the positive binding poten- 
tial. For circular orbits the integrals of motion can be written 
in function of the radius r. 



r Alb 
2 ar 



dr 



(9) 



(10) 



which can be solved to obtain EdL) and LdE). All dynamical 
properties can be derived from the DF, such as the anisotropic 
velocity moments 



dv,. 



F(E,L)vfvl'"^^ dvr, 



(11) 

with Mtot the total mass. If the system is self-consistent, then 
the density p{r) = /ioo('") is connected to the potential via the 
Poisson equation 



1 d / , dV' , 



(12) 



Furthermore, the second-order moments determine the ra- 
dial and tangential velocity dispersions /i2o('") = Pf^(r) and 
IJ-oiif) = 2pal(r), and the velocity anisotropy profile 



/3(r)=l- 



(13) 
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A DF is only ful ly determine d when all velocity moments dTTT i 
are known (De jonghelll987h . However, it is akeady far from 
trivial to obtain any non-negative DF that adequately gener- 
ates a given density p(r) and anisotropy profile (3(r). In this 
paper, we will focus our attention to one specific family of 
halos, namely the Dehnen-McLaughlin Jeans models, and we 
demonstrate that a family of DFs can indeed be constructed 
for these systems. 

2.2. The Dehnen-McLaughlin halos 

In the c ontext of theoretical dark matter studies, the model 
derived by lDehnen & McLaug hlin (2005) is of particular in- 
terest. We summarize their main results in this section. In- 
stead of fitting a parameterized density profile to A^-body sim- 
ulations, they investigated the solution space of the Jeans 
equation to search for models that explicitly obey the power- 
law behavior of the pseudo phase-space density. With the ex- 
tra condition of a linear f3-j relation they found a critical 
solution that satisfies the condition 



(14) 



with Ts a scale radius. In the remainder of this paper, we adopt 
the common value e = 3 and use the notation Qrir) = p/crl- 
Dehnen & McLaughlin derived for this case the exponent 



phase-space density, two scaling constants i.e. the total mass 
Mtot and the scale radius r^, and the asymptotic anisotropy 
parameters /3o and j3oo ■ The authors also noted the remarkable 
property that the shape of the density profile (and hence the 
gravitational potential) only depends on /3o and not on /3oo ■ 

While the Dehnen-McLaughlin Jeans models are derived 
from theoretical considerations, they also fit adequately 
galaxy-sized and c ertain cluste r-sized halos generated by A^- 
body simulations (iDiemand et al. 2005; Mer ritt et all (2006). 
But as mentioned in the previous Section, the density and ve- 
locity dispersions alone do not determine the complete dy- 
namical state of dark matter systems. We therefore aim to 
incorporate these profiles into self-consistent dynamical mod- 
els, described by non-negative DFs. 

3. THE MODELING TECHNIQUE 

To find DFs that describe the Dehnen-McLaughlin Jeans 
models, we adopt the mathematical framework that we de- 
rived in Paper I. These tools enabled us to construct a fam- 
ily of components with the anisotropy profile Q, by means 
of the powerful augmented density concept. We summarize 
these results, and we demonstrate how to build a linear com- 
bination of these components using a quadratic programming 
(QP) technique (Dejonghe 1989). In this manner, we can fit a 
dynamical model to a given halo. 



77 = 



= 77+-, 
' 2' 

4-2A) 



(15) 
(16) 



with (3q the central velocity anisotropy. The corresponding 
potential-density pair is fully analytical, given by 



-0(r)= B_ 



1 i-A)^i 



(17) 



4 + 77-2A, M^^_,„ ^^,^-(,^-,0)/, (^g^ 

where x = r/r^, By{a,b) is the incomplete beta function and 

7+10/3o 



70= ■ 



31-2/3o 



(19) 



(20) 



The density can be equivalently written in terms of the 
slope 7(r), which has the same elegant form as the velocity 
anisotropy profile (3(r) 



7(r) = 
/3(r) = 



(21) 
(22) 



l+x*) 

Finally, the authors derived the corresponding velocity disper- 
sions 



a?(r) = 



1 



4 + r;-2/3oo r, 
1 



-X 



l+x'' 



crj(r) = al(r)=-aUr)={l-m) <K{r). 



,(23) 



(24) 



Eqs. (fT4b-(l24l) characterize the Dehnen-McLaughlin models; 
we refer to their paper for more details. These systems are 
determined by five parameters: the exponent e in the pseudo 



3.1. The augmented density concept 

Since our dynamical models have to reproduce the mo- 
ments ( fTSl l. ( |23] ) and ( |24] |. we first seek DFs that are specifi- 
cally designed for this task. This c an be do ne by introducing 
the augmented densities piip^r) CDeionghe 1986), which ex- 
tend the densities to explicit functions of both the radius and 
the gravitational potential. Like the DF, the augmented den- 
sity uniquely determines the dynamical state of a spherical 
equilibrium system. Both functions are connected by the re- 
lation 



p(7/;,r) = 2^Mto, / AE 
Jo Jo 



2(i,-E) 



F(E,rvT) 



= dv 



T- 



(25) 



If we define that all functions are zero when their arguments 
lie outside their physical bounds, we can use the Laplace- 
MelUn transforms 



/>+oo 

C M {F}= e'^^dE / L^-^F{E,L)AL 

s^C L^A Jo Jo 

C M {~p}= I 



(26) 



^+00 p-\-oo 



d^ / r^-'p{iP,r)Ar, (27) 




with the inverse transforms 



_2 /-Co+'oo 



F(E,L)=-^ e«^d£/ L-^ C M{F}dL,(28) 

47r^./f„_;_ JAo-!00 L^X 

e«'^dV'/ C M{p}dr. (29) 



Combining these expressions with (1251 1. we obtain 



£ M {p} = 27rMtot / 

I'^i Jo 



dTA 



X dE 
Jo Jo 



F{E,L)dv} 
^2(^-E)-vj 



(30) 
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Interchanging the integrals yields the formal relation 

2A/2 ^(3-A)/2 



E^i Z^A M,„,(2^)3/2 r{l-^)Ai 



£ Mjp}. (31) 



The main strength of the augmented density formalism is the 
abiUty to impose very specific conditions on these functions, 
to obtai n our objective. More precisely, we demonstrated in 
iPaper ll that separable functions of the form 

p(iP,r) = f(iP)g(r), (32) 

are particularly interesting, since in this case g(r) is directly 
related to the velocity anisotropy profile. 



/?('") = -:^ -7; — ('")• 



(33) 



2dlnr 

With this formalism it is therefore at least formally possible to 
generate a DF for a dynamical system with a given potential 
ipir), density p(r), and velocity anisotropy profile /3{r). If we 
postulate an anisotropy profile of the form (|5]), the modeling 
procedure reduces to the computation of f(4>) such that 



p(r) = f{ij(r)) 
with < (5 1 and 



-2/3o 



^25 \ 



.2(5 



f3s = 



f3o-f3. 



(34) 



(35) 



Evidently, the exact solution /(ip) is generally of numerical 
form, except for a few special cases. This creates a prob- 
lem to recover the DF from the augmented density, because a 
direct inversion of the Laplace-Me llin transformat ion in OTb 
is in general numerically unstable (Dejonghd ll986) . Indeed, 
the double integration in ( |25] | smooths features in F(E,L) and 
the inversion procedure of determining F{E,L) from pi'ip,r) 
has the delicate job of unsmoothing the information contained 
in p(^p, r). Consequently, a direct inversion can only be per- 
formed safely for sufficiently simple forms of f(^). 

We therefore take an alternative approach. Instead of 
solving ( I34I 1 directly, we will approximate the given den- 
sity profile by a linear combination of simple base functions 
Pi(4>ir),r) for which the Laplace-Mellin inversions are ana- 
lytical. The corresponding DF is then simply the same linear 
combination of the associated base DFs Fi(E,L), resulting in 
an analytically tractable function. This can be achieved by 
a least squares fit to a set of density data points, defining a 
quadratic programming problem in the unknown coefficients. 
Various au thors have successfully used a similar modeling 
techii i aue fcuiiken & Merrifieldlll99"l iMerrifield & KuijkenI 
Il994t iGerhard et al.lll998l) . AOP-algorithm. suited for this 
task, has been developed in our department ( Dejonghe. 1989.) , 
which we outline in the following Section. 

3.2. The quadratic programming procedure 

Consider a given potential ipir), an anisotropy profile /3(r) 
of the form (|5]) and a set of M density data points Pobsirm), 
m = 1, . . . ,M. To model these data, we first construct a library 
of A^iib base functions 



p,#,'-) = .//(V') 



-2/3o 



^25 \ 



.2S 



(36) 



with fii^) sufficiently simple to compute the associated dis- 
tributions Fi(E,L). From these components, we can extract 
the corresponding values 



Pi(rm) = Pi{ip(rm),i 



m = 1, 



(37) 



Our aim is now to construct a linear combination of com- 
ponents from this library that provides an adequate fit to the 
given data. Such a fit can be obtained in a statistically mean- 
ingful way by minimizing the quantity 



1 ^ 

m=l 



N 

1=1 



Pobs('"m)- > aN,iPi(r„) 



(38) 



which is a quadratic function of the coefficients a^,,, con- 
sequently defining a quadratic programming problem. The 
data points are given equal weights by setting the constants 

To find such a set, we use an iterative algorithm (Uejonghg 
Il989l) . With this method, a set of components is successively 
built in steps. In the first step, Nnh x^-minimizations are 
performed using in turn each component of the library. The 
base function that yields the lowest x^-value (hereafter de- 
noted x^) is retained as the first element of our best-fitting 
set, with corresponding coefficient a 1.1. In each next itera- 
tion, the functions in this set are preserved, while the coeffi- 
cients are allowed to vary. The set is subsequently extended 
by adding the component from the library that yields the most 
improvement of the fit, minimizing over all coefficients. In 
other words, suppose we have obtained the best-fitting set 
of A^- 1 base functions, with indices ni,... ,niv-i. Then, we 
add in turn the remaining Nnh -N+1 components from the 
library, and calculate the coefficients for each combination 
by means of the -minimizations 



X(n«) " 



^ M / N \ ' 

= min — y^w„, \ Pobsirm) -y^aip„.(rm)\ 

a\,....,an M — ' \ ^ — ' / 
m=I \ 1=1 / 

with «A, e {1, . . . ,A^iib}\{«i, . . . ,«iv-i}. 



From these N\ib-N+ 1 values, we determine the best fit 



(39) 



(40) 



and add the base function with index Mmin to the best-fitting 
set, denoting n^i = By renaming the indices of this 

set, we thus obtain a linear combination of A^ base functions 
Pi(^,r) with coefficients a^j and a goodness of fit Xn- The 
corresponding DF is then simply 



FiE,L) = Y,aN.iFiiE,L). 



(41) 



The QP-algorithm allows additional linear constraints on the 
coefficients. In particular, we impose upper and lower bound- 
aries 

max; 

VA^; /= l,...,A^. (42) 

These constraints are not necessary, but they greatly reduce 
the computational cost in the calculation of the DF. The reason 
for this is straightforward: if the components are computed 
with numerical errors (5,7^ then the total numerical error of the 
DFis 



(5F(£,L): 



S,=i \a\N.iSiF(E,L) 



(43) 



The higher the absolute values of the coefficients |fl|iv,,, the 
smaller the errors SiF, need to be to obtain a given SF, which 
increases the computational time. Sensible boundary values 
( |42] | enable efficient calculations of the DF, while maintaining 
satisfactory fits. 
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The resulting DF also needs to be physical, i.e. non- 
negative everywhere in phase space. We found that all our 
DFs automatically satisfy this condition without imposing ex- 
plicit constraints. 

This procedure has several advantages. The resulting DF 
remains analytically tractable, which also simplifies the com- 
putation of all the moments of this dynamical model. Further- 
more, only a limited number of data points are required, rather 
than the entire density profile. The same algorithm can also 
be applied to data extracted from simulations or observations, 
and a variety of different moments besides the density can be 
used in the fitting procedure. 

Finally, since all components have a priori the desired /3(r), 
their linear combination will automatically generate the same 
anisotropy profile. This is a very significant benefit. Indeed, 
practice has shown that it is particularly difficult to construct 
radially anisotropic systems from b ase functions t hat are too 
simple (such as Fricke components, iFrickd (11952 !)) with dif- 
ferent constant anisotropics. In that case an additional fitting 
is required to the second-order moments. But more impor- 
tantly, because radial orbits influence the density both at small 
and large radii, the summation of components with different 
anisotropics results into a delicate fine-tuning to obtain ad- 
equate fits to both the density and the velocity anisotropy. 
These problems are avoided if the components already have 
the correct anisotropy profile, and although it is more intri- 
cate to design such functions, this approach greatly reduces 
the complexity of the quadratic programming procedure. 

3.3. The base functions 

We are now left with the construction of a library of ade- 
quate base functions r). The success of our modeling 
is largely determined by this library. As stated above, the 
components have to generate the anisotropy profile (|5]), while 
being sufficiently simple to retrieve the corresponding DFs 
Fi{E,L) from eq. dSTl i. In addition, the subsequent densities 
p,(rm) need to be able to reproduce the specific characteristics 
of the given data to obtain a satisfactory fit. In the particular 
case of the Dehnen-McLaughlin halos ( fTsl l. this implies that 
the components should incorporate the central density cusp 
and the asymptotic density slope at large radii, using the given 
potential ( fTTI ). In [Paper H we derived a family of base func- 
tions that meet all these requirements. Consider the family of 
augmented densities 



-2/3o 



1 + - 



.25 



where poi are normalization constants such that 



47r 



Pi(-4)(r),r)p-Ar= 1, 



(44) 



(45) 



and denotes the depth of the (finite) potential well, ?/)o = 
V'(O). The normalization constants are chosen such that the 
sum of the coefficients of a good fit should approximate the 

total mass of the input data, i.e. X^ti — ■'^tot- 
Apart from the four fixed parameters that determine the 
anisotropy profile (l22T i. the /(V')-part of these functions con- 
tains three additional free parameters pi, qi and i,, which re- 
spectively determine the asymptotic behavior at infinity, the 
inner slope, and the transition rate between these two regions. 
They satisfy the conditions p, + 2/3oo > 3, qi ^ and > 0. 

All other dynamical properties can be calculated from these 
augmented densities. The augmented higher-order moments 



are derived in Appe ndix lAl and we demonstrated in the Ap- 
pendices of iPaper j that for these /9,(V', r) the inverse Laplace- 
Mellin transforms can be performed analytically in eq. (|3TI ). 
The corresponding distribution functions can be expressed as 
a series of Fox //-functions (lFoxlll961l) : 



M, 



tot 



X^(-l) 



qi\ r(i +/>/+;«,■) 



5T{-I3s) 



1.1 

2.2 



L2 



2r?£ 



[ 5' 

which can be written as a double series 



13^ i 
5 ' S 



,(0,1) 



,(46) 



POi 



m+pi+jsi) 



,.0 s - xT (;,,■ + ;•.,- i + A) r (1-/3,) 
where we used the auxiliary notation 

/3o-M for L^<2rlE 



2rlE 



for U > IriE. 



(47) 



(48) 



The double summation J^j J2k '-^'^ computed by chang- 
ing the indices to J^i 12j+k=i^ ^^^^ '■^^ inner summation be- 
comes a finite sum of / + 1 terms for each value of /; the index 
I is increased until the total sum alters by less than a required 
numerical error SiF,. Due to the double summation, the com- 
putational time i s an inve rse quadratic function of 5,/] . 

We proved in iPaper II that these base DFs are continuous 
and non-negative everywhere in physical phase space. More- 
over, we showed that they gene rate exactly the generalized 
Plummer potential-density pairs (I VeltmannI 1 9791) 

^gpW=Ti^^, (49) 



Pgp('") = 



(r? + /')'/" ^ 
(l + 77)M tot 
4n 



V 



(50) 



with a linear - 7 relation. Since these systems are closely 
related to the Dehnen-McLaughlin halos, this is a promising 
result for the success of the quadratic programming routine. 

The Fortran source code with the DF base functions and 
the augmented moments is available on request. 

3.4. The library of components 

Every given Dehnen-McLaughlin halo requires a specific 
component library. In particular, the parameters qi are 
constrained by the potential. If we examine the asymptotic 
behavior of the Dehnen-McLaughlin potential (fTTb in more 
detail, we find 



V'(r) - i/;o-fl/"-"''^°>/'' + --- for r^O, 

ip{r) ^ — for r^oo. 

r 



(51) 



Introducing these asymptotic expansions in the expres- 
sion (l44l) we find for the inner and outer slopes of the density 

^-2A+?,(ll-10/3o)/9 for ^^0, 



PiWr),r) 
Pi{ip{r),r) 



for r 
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Fig. 1 . — The obtained xjj for th^ six QP-models, explicitly as a function 
of the number of components in the fit. The different curves correspond to 
/3oo = 0, 0.2, 0.4, 0.6, 0.8 and 1, with a grayscale ranging from black ($00 = 0) 
to lightgray (/3oo = 1). 

Evidently, the parameters pi stipulate the density slope at 
large radii. Because the models fall as r~''°^, no components 
can be used in the fitting routine that fall less rapidly. Using 
eq. (l20l i. this puts a boundary on the p,, 

31-2/3o-18/3oo _ 



Pi 



(/?0,/?oo). 



(53) 



Conversely, the density slope at small radii depends on the 
parameters The density diverges toward the center as r"^°, 
and we cannot use components in the fitting routine that have 
a steeper slope. Thus we obtain from eq. ( |20] | 



qi > 



7-8/?o 

ii-iOA) 



^ ?min(/3o)- 



(54) 



So if a fit to a halo has at least one component with parameter 
Pmin and one with qmin, this fit has the same slope as the given 
density at small and large radii. 

Finally, the parameters i, have a similar role as 6, in the 
sense that they control the transition rate between the inner 
and outer density slopes. Their value can be chosen freely, 
but we found that excellent results are obtained with a single 
fixed value 

Si = 25 = r], (55) 

for all components. This is th e same choice as the generalized 
Plummer models in 'Paper f. It also simplifies the computa- 
tion of the DFs (|42ll-(I^, and it further facilitates the fitting 
process, leaving only pi and qj as free parameters. 

4. RESULTS 
4.1. The minimization 

Now that we have derived the necessary mathematical tools, 
we can present the results for the Dehnen-McLaughlin halos. 
Without loss of generality, we can work in dimensionless units 
G = M(ot = '"s = '"a = 1^ and we limit ourselves to e = 3. Con- 
sequently, the models are determined by the anisotropy pa- 
rameters /3o and (ioo ■ Although we are able to generate mod- 
els with arbitrary values for these parameters, realistic dark 
matter halos are nearly isotopic near the center and radially 
anisotropic at large radii, so that we concentrate on six repre- 
sentative models with /3o = and /3oo = 0, 0.2, 0.4, 0.6, 0.8, 1. 
We verified that the modeling procedure works equally well 
for models with non-zero values of Po. Finally, it is evident 
from eq. (fTSI l that /3o = sets the parameters si = 25 = 77 = 4/9. 

As we demonstrated above, the very specific form of the 
base functions (l44l i simplifies our QP-algorithm considerably 
for these models. Only the parameters pi and qi remain to 
construct a library of components, and we have found that 



icr* 
icr' 

ID-" 

-10-'° 

-Iff' 

-icr* 

-lET^ 



1 1 n-PT 1 ■ ■ ■ 1 ■ 


' ■ ' 1 ' ■ ' ' 1 ■ ' " 


T ■ ' ■ ' 1 ■ ■ ' ' 


, ^^^J , . 


"""^^ 

, . 


-1 ■ ■ 










■ ■ ■ ■ ■ ■ ■ ■ ° I ■ u 





lD-= 



Iff' 



rjri 



ID* 



Fig. 2. — The 10 individual components of the fitted density for the QP- 
model with /3o = and Poo = 0.4. Their sum is the QP-density (black thick 
curve), fitting the 25 data points (dots). 



only 30 components are sufficient to yield excellent fits. The 
parameters p, take five values, ranging from Pmmif^iPoo) to 
10 or 12, depending on the model. The parameters qi take six 
values from <7min(0) to 0. We list these values for each model 
in the headers of Table [T] 

We extract M = 25 values of the density dTsl l at radii r,„, dis- 
tributed logarithmically between IQT^ and 10** r,, that serve 
as input data Pohsifm) in each QP-procedure. Evidently, this 
range is much larger than the virialized region in A^-body sim- 
ulations. This larger range is therefore not intended to be real- 
istic, but rather to demonstrate that our models are accurate to 
arbitrary distances. Furthermore, this makes it possible to cre- 
ate discrete equilibrium systems from the DFs, by means of 
Monte Carlo simulators, that trace very closely the Dehnen- 
McLaughlin halos. After calculating the densities ( l44l i of ev- 
ery library component at these radii pi{ip{r,„), r,„), we can per- 
form the QP-algorithm for the six values of f3oo, constructing 
iteratively the best-fitting linear combination ( l38b of com- 
ponents with additional constraints of the form ( l42b . 



100 s$ gnj < 100, 



VA^; /= 1,. 



(56) 



We show the results for the six models in Table [T] The 
columns list the components of the fits. Every xli denotes the 
goodness of fit of the best linear combination of the compo- 
nents in columns 1 to A^. Each component is determined by 
the parameters pi, and which in turn define the normal- 
ization constants po/ from eq. ( l45b . The coefficients are only 
given for the final fit with 10 components, i.e. oio ,. It can be 
checked that for each model 



10 

E 

/=i 



am.i —Mtot = 1, 



(57) 



indicating excellent fits. Combining this result with eqs. ( l43T l 
and ( |56] |. it can be seen that if A^ = 10, the numerical errors 
of the base functions (5,/v need at most be a factor lO^* smaller 
than a given error SF, allowing efficient computations of the 
DF with sufficient accuracy. 

The resulting xli for each model are also displayed in Fig.[T] 
Evidently, A' = 10 components are more than sufficient to ob- 
tain very accurate dynamical models. As an example, Fig.|2] 
shows the 10 individual components of the QP-model with 
(3oo = 0.4. Although this fit has the highest Xio of our set, its 
total density is a very close approximation to the given data 
over the entire range in radius. 
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TABLE 1 

Components of the six QP-models. 





1 


2 


3 




4 


5 


6 


7 


8 


9 


10 




/3o = 


= 0.0, 13^ -- 


= 0.0, s = 


4/9, 


<5 = 2/9 


Pi = 


3.44, 4, 5, 7.5, 10 




-0.63, -0.5, -0.4, 


-0.3, -0.15,0 




"lOJ 


0.2047 


-0.0380 


-0.5411 




1.6652 


0.0414 


-3.8806 


0.1312 


2.5207 


-0.2684 


1.1652 


POi 


0.0053 


1.9783 


0.0009 




2.5540 


0.4890 


3.0677 


0.0384 


3.6728 


6.1920 


0.0008 


Pi 


4.0000 


10.0000 


3.4444 




10.0000 


7.5000 


10.0000 


5.0000 


10.0000 


10.0000 


3.4444 


qi 


-0.6364 


-0.6364 


0.0000 




-0.5000 


-0.6364 


-0.4000 


-0.6364 


-0.3000 


0.0000 


-0.6364 




0.44- 10" 


0.17 • 10" 


0.67 ■ 10"' 




0.33 - 10"' 


0.24- 10-2 


0.17 - 10"2 


0.53 - 10"^ 


0.36- 10^ 


0.48 - 10"5 


0.82-10-' 




/3o = 


= 0.0, /3oo = 


= 0.2, s = 


4/9, 


5 = 2/9 


Pi = 


3.04, 4, 5, 7.5, 10 


li = 


-0.63, -0.5, -0.4, 


-0.3, -0.15,0 






-0.1534 


15.5467 


-0.0006 




0.4549 


-0.2707 


74.2329 


-88.6417 


0.5880 


-0.7724 


0.0165 


POi 


0.0238 


0.0017 


3.4087 




4.4609 


11.3269 


0.0014 


0.0015 


8.6356 


5.4102 


0.9797 


Pi 


4.0000 


3.0444 


10.0000 




10.0000 


10.0000 


3.0444 


3.0444 


10.0000 


10.0000 


7.5000 


1i 


-0.6364 


0.0000 


-0.6364 




-0.5000 


0.0000 


-0.6364 


-0.5000 


-0.1500 


-0.4000 


-0.6364 


AN 


0.48- lO" 


0.27 ■ lO" 


0.43 ■ 10"' 




0.93 ■ 10"- 


0.13 - 10"2 


0.88 - 10"^ 


0.10- 10"^ 


0.95-10-5 


0.48 - 10"5 


0.32- 10"'' 






= 0.0, /3oo 


= 0.4, s -- 


= 4/9, .5 = 2/9 


Pi - 


= 2.64, 4, 5, 8, 12 


qi = - 


-0.63, -0.5, -0.4, ■ 


■0.3,-0.15,0 




Olo.i 


8.2375 


-1.4318 


-0.0467 




0.2297 


2.2705 


-0.0111 


0.1299 


10.3097 


-0.3108 


-18.3771 


POi 


0.0832 


0.0032 


10.9019 




14.9246 


0.0028 


2.4301 


23.3335 


0.1015 


18.6978 


0.0934 


Pi 


4.0000 


2.6444 


12.0000 




12.0000 


2.6444 


8.0000 


12.0000 


4.0000 


12.0000 


4.0000 




-0.6364 


0.0000 


-0.6364 




-0.5000 


-0.6364 


-0.6364 


-0.3000 


-0.4000 


-0.4000 


-0.5000 


AN 


0.48 ■ 10" 


0.18-10" 


0.28 ■ 10-' 




0.93 - 10"2 


0.42 - 10"2 


0.14-10"^ 


0.36- 10~* 


0.13 - 10"* 


0.68 - 10"5 


0.24 - 10"* 




/3o = 


= 0.0, /3oo = 


= 0.6, s = 


4/9, 


<5 = 2/9 


Pi = 


2.24, 3.5, 5, 8, 12 


qi = 


-0.63, -0.5, -0.4, 


-0.3, -0.15,0 




aio.i 


0.3552 


2.0738 


-0.0094 




0.0086 


16.3498 


-17.4861 


0.0731 


-0.3019 


-0.2109 


0.1477 


POi 


0.1222 


0.0063 


16.7922 




23.2672 


0.0053 


0.0055 


4.3156 


0.1372 


8.0592 


10.5109 


Pi 


3.5000 


2.2444 


12.0000 




12.0000 


2.2444 


2.2444 


8.0000 


3.5000 


8.0000 


8.0000 


li 


-0.6364 


0.0000 


-0.6364 




-0.5000 


-0.6364 


-0.5000 


-0.6364 


-0.5000 


-0.3000 


-0.1500 


■xl 


0.44 • 10" 


0.15 ■ lO" 


0.18-10"' 




0.67-10"- 


0.31 - 10"^ 


0.16-10^ 


0.13 - 10"* 


0.10- 10"* 


0.17 - lO""" 


0.20- 10"* 




ft 


= 0.0, /3oo 


= 0.8, S-- 


= 4/9, S = 2/9 


Pi - 


= 1.84, 3, 5, 8, 12 


li = - 


-0.63, -0.5, -0.4, ■ 


■0.3,-0.15,0 




aio.i 


-0.0321 


-6.1421 


0.0008 




-0.0013 


24.5830 


-100.0000 


-0.0163 


-0.0926 


40.5233 


42.1774 


POi 


0.1791 


0.0121 


25.3820 




35.5933 


0.0101 


0.0109 


2.3877 


0.2367 


0.0105 


0.0116 


Pi 


3.0000 


1.8444 


12.0000 




12.0000 


1.8444 


1.8444 


5.0000 


3.0000 


1.8444 


1.8444 


1i 


-0.6364 


0.0000 


-0.6364 




-0.5000 


-0.6364 


-0.4000 


-0.3000 


-0.3000 


-0.5000 


-0.1500 


Xn 


0.40 • 10" 


0.11 ■ 10" 


0.11 - 10"' 




0.46- 10"2 


0.22 - 10"2 


0.66 - 10"5 


0.26 - 10"5 


0.13- 10"5 


0.47 ■ 10"* 


0.24-10"* 




ft 


= 0.0, /3oo 


= 1.0, S-- 


= 4/9 


5 = 2/9 


Pi-- 


= 1.44,3, 5, 8, 12 


* = - 


-0.63, -0.5, -0.4, ■ 


■0.3,-0.15,0 




aioj 


-5.2047 


0.1086 


62.6962 




40.1800 


-0.0002 


0.0003 


-96.6969 


-0.1234 


0.0273 


0.0128 


POi 


0.0192 


2.9967 


0.0201 




0.0213 


37.6938 


44.2497 


0.0207 


3.7237 


7.8118 


0.8543 


Pi 


1.4444 


5.0000 


1.4444 




1.4444 


12.0000 


8.0000 


1.4444 


5.0000 


5.0000 


3.0000 


1i 


-0.6364 


-0.6364 


-0.5000 




-0.3000 


-0.6364 


0.0000 


-0.4000 


-0.5000 


0.0000 


-0.1500 


xl 


0.35 ■ 10" 


0.36 - 10-' 


0.10-10"' 




0.24- 


0.41 - 10"^ 


0.16-10"^ 


0.13-10^ 


0.45 - 10"^ 


0.89-10"^ 


0.22 - 10"** 



Note. — Our six models are determined by the anisotropy parameters ft = and /3ac- With the QP-algorithm, we built linear combinations up to 10 
components, each chai'acterized by s, 5, p;, and po,. The parameters s and 5 are the same for each component, while p, and qj are selected from a library of 
30 components. The five p; values and six qj values in each Hbrary are given in the headers. Every A'th column lists the value xl of 'he best fit with the first N 
components. The coefficients are given for the final fit with 10 components, i.e. «io,i. 



4.2. The moments 

Fig.[3]displays several moments for our six models, with 10 
components. The top row shows the density p(r), the pseudo 
phase-space density QAf), and the /3- 7 relation. Below each 
graph, we calculated the residual errors between the QP-fits 
and the theoretical curves, i.e. for each profile /(r) we have 

/obs('-)-/qp('") 



A/(r): 



./obs('") 



(58) 



As can be seen, the relative errors on the densities are less 
than 10"^ along 7 orders of magnitude in radius, and the ex- 
act asymptotic slopes of the models ensure excellent fits even 
beyond this range. The power-law trend of QAf) and the 
/3 - 7 relations are also reproduced very accurately with er- 



rors ^ 10"^. Note that the small offset between the pseudo 
phase-space density profiles for the different models is due to 
the dependence of cr, (r) on (3oo- 

In the central row, we display the velocity dispersion pro- 
files ar{r), crg{r), and the anisotropics (3(r). It is striking 
that, while the dispersion data were not used in the fit, the 
deviations of these moments from the theoretical values are 
even smaller, less than 5 - 10"^. Evidently, since the models 
have the anisotropics (|22] | by construction, the (3(r) profiles 
are exact, without errors. Note also that all tangential ve- 
locity dispersion profiles aeir) intersect at a common radius 

^=^,(9/11)'/". 

While the density and dispersions are defined by the 
Dehnen-McLaughlin halos, the higher-order moments are de- 
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Qr{r) 
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itf ifl* 
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m' id-'' 1^ 1^ ifl* IB' is' IB''' 1^ 1^ 1^ 10* 
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Fig. 3. — The most important moments for our set of representative models with 10 components. Top row: the density p(r), the pseudo phase-space density 
Qr(r) and the /3-7 relation. Below each graph, the relative errors with respect to the theoretical profiles is shown. Middle row: the radial velocity dispersion 
(T,-(r), the tangential velocity dispersion (Tg(r), and the anisotropy /3(r), also with the relative eiTors. Bottom row: the radial kurtosis Kr(r), the tangential kurtosis 
Kg{r), and fourth-order anisotropy /34(r). The models and grayscaling are the same as in Fig.[T] 

termined by the QP-models. The fourth-order moments (see 
Appendix |A]i allow us to derive the radial and tangential kur- 
tosis. 



«e(r)=^(r)-3, 

Co 



and the fourth-order anisotropy 



/34W=l-|4f(r). 



(59) 
(60) 

(61) 



Interestingly, as a result of the separable form of the aug- 
mented densities, we find that the /34(r) profiles are only a 
function of the /3(r), 

fiAir) = (3 - m) + (A) " m) {Poo - fi{r)) . (62) 

These profiles are shown in the bottom row of Fig. |3] The 
kurtosis values describe the non-Gaussianity of the veloc- 
ity distributions at a certain radius. Our radial kurtosis val- 
ues are very large in the center, which indicates that the 
Vf-distributions are very peaked (leptokurtic) at small radii. 
The K,(r) curves decrease rapidly as a function of radius: 
they reach zero at radii between 0.26-0.36 and become neg- 
ative at larger radii, leading to flat-topped (platykurtic) radial 
velocity distributions. This behavior is in accordance with 



A^-bo dy simulations (iKazantzidis et al.l l2004t IWojtak et al.l 
l2005h . Clearly, the value of /3oo has little influence on the ra- 
dial kurtosis, as in the case of Qr{r). In contrast, the tangential 
kurtosis no{r) curves do depend significantly on Poo- All vg- 
distributions are highly peaked at small radii. For j3oo < 0.4 
the tangential kurtosis decreases to slightly negative values, 
i.e. at larger radii the tangential velocity distributions become 
slightly flat-topped. Models with Poo ~ 0.4 have nearly Gaus- 
sian V0-distributions for r > r^. If Poc > 0.4, the Ke{r) profiles 
reach a minimum value and increase again for larger radii. 

4.3. The distribution functions 

The six top panels of Fig. |4] show the distribution functions 
F(E,L) of our radially anisotropic systems with /3o = and 
Poo = 0, . . . ; 1, expressed as logarithmic isoprobability con- 
tours and a logarithmic color gradient in the integral space, 
with L scaled to L^^iE), denoting the angular momentum of a 
circular orbit with energy E. All models are clearly physical, 
i.e. the DFs are non-negative everywhere. This means that 
the Dehnen-McLaughlin Jeans models can actually be real- 
ized by full dynamical models. Moreover, contrary to the 
Osipkov-Merritt models, these functions fill the entire inte- 
gral space. In the isotropic case, the contours are horizontal 
(no dependence on angular momentum), and their orientation 
alters gradually with increasing /3oo in an intuitive way, as 
orbits with high eccentricities (i.e. low angular momentum) 
become more abundant. 

The DFs describe the probability distributions of particles 
in phase space, but not in the integral space. It is more physi- 
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Fig. 4. — The phase-space distribution functions for our set of models with 10 components. The six top panels show F{E,L) as isoprobability contours in the 
integral space. The energy is scaled to the central potential and the angular momentum is scaled to the angular momentum Lc (£") of a circular orbit with energy 
E. The contour levels and the coloring are scaled logarithmically. The six bottom panels show the orbital distribution functions N [E ,L/Lc(E)) for the same 
models. 
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cally meaningful to consider the true orbital distributions 

N {E,L/L,(E)) = F(E,L)g(E,L)L,(E), (63) 
with the so-called "density of states" function 

dr 



giE,L)=16TT'L 



r_ ^2(i,(r)-E)-Lyr 



(64) 



where r_ and r+ are the pericenter and apocenter of an orbit 
with energy E and angular momentum L. In other words, the 
functions A^(£',L/Le(£')) express the likelihood of an orbit 
with energy E and scaled angular momentum L/LdE). 

The results are displayed in the bottom panels of Fig. |4]as 
logarithmic isoprobability contours and a logarithmic color 
gradient in the integral space. For high binding energies, all 
orbital distributions contain increasing probabilities toward 
circular orbits (high angular momentum). For low binding en- 
ergies, i.e. at large radii, the abundance of orbits with low an- 
gular momentum gradually increases from the isotropic case 
to models with high values of Poo- 

4.4. The marginal distributions 

We conclude the discussion of our Dehnen-McLaughlin 
DFs with an analysis of the marginal distributions. The dif- 
ferential energy and angular momentum distributions are the 
integrals of the orbital distiibutions. 



f.L,(E) 

N{E)= I N(E,L)dL, 
NiL)= I NiE,L)dE, 



with 



N(E,L): 



L,iE) 



N{E,L/L,(E)). 



(65) 
(66) 

(67) 



These curves are displayed in Fig. |5] The differential energy 
distributions are all monotonously decreasing functions of E. 
It is striking that these profiles are almost identical, regardless 
of the anisot ropy 0^^ . Thi s result reinforces previous dynam- 
ical studies (iBinnevll 19821) . and suggests that p(r), Qr(r) and 
Krir) are linked to a universal differential energy distribution, 
independent of (3oo, caused by the same physical processes. 

In contrast, the angular momentum distributions depend on 
/3oo, most notably for high radial anisotropics. For increas- 
ing values of Poo, t he fraction of o rbits with low angular mo- 
mentum increases. ' Bullock et alJ fcOOli) proposed a univer- 
sal form for the integrated angular mo mentum distribution in 
dark m atter halos M(L). Alternatively, ISharma (feSteinmetzl 
(l2005h found a differential distribution 



1 



L°r(fl) 



-L^-'e"^/-^". 



(68) 



Our models indicate a similar profile, with a > 0.9, although 
the functions ( |68T l fall steeper than ours as L increases. 

5. CONCLUSIONS 

In this paper we presented a set of anisotropic distribution 
functions F(E,L) for the Dehnen-McLaughlin Jeans models. 
We constructed these DFs as a linear combination of base 
functions, which we fitted to the halo density by means of 



a quadratic programming algorithm. The base functions were 
specifically designed for this task, derived from a separable 
augmented density that generates exactly the general four- 
parameter velocity anisotropy profiles (|5]). We demonstrated 
that the resulting fits are very accurate, from the center to large 
radii, far beyond the range of A^-body simulations. 

This method has several advantages. The DFs can be writ- 
ten as a sum of a double series, without numerical integrations 
or inversions. Consequently, the DFs and all subsequent mo- 
ments can be computed with high accuracy. Moreover, the ad- 
vanced form of the base functions makes the QP-fitting very 
fast, requiring only a small number of library components. 

In this manner, we have constructed a family of dynami- 
cal models that incorporate the observed features in A^-body 
simulations. We summarize their properties: 

1 . The models a priori have the universal properties en- 
countered in A^-body studies of dark matter halos, 
w hich formed the building bricks of the Jeans models 
bv'D ehnen & McLaughhnl(l2005h : they generate a uni- 
versal density profile, a power-law pseudo phase-space 
density QAf), four-parameter velocity anisotropy pro- 
files with arbitrary values of Pq and /3oo, and a linear 
/3 - 7 relation. In particular, we analyzed six models 
that are isotropic in the center and radially anisotropic 
at large radii. 

2. The DFs are physical, i.e. they are non-negative over 
the entire phase space. This means that the Dehnen- 
McLaughlin Jeans models can actually be realized by 
self-consistent dynamical systems. In addition, the DFs 
are continuous and smooth functions and, contrary to 
the popular Osipkov-Merritt models, they fill the entire 
accessible phase space. 

3. The energy distributions are monotonously decreasing 
functions, and like the radial kurtosis profiles, these 
functions are nearly independent of /3oo- This suggests 
that p(r), Qrir), Krir), and A^(£') have a universal form, 
caused by the same physical processes. 

In addition, we have also been able to generate dynamical 
models with a non-linear /3-7 relation, i.e. S y^i]/! and ra ^ rj, 
albeit with higher xjj values. From these successful results for 
the Dehnen-McLaughlin halos, we expect equally adequate 
fits for other Zhao models, such as the NFW halos. Other pro- 
files, such as the Sersic-type densities, might require a modi- 
fied family of base functions. Such a systematic study could 
unravel more hitherto hidden properties of dark matter halos, 
and the various connections between these characteristics. 

Our set of dynamical models is not only useful for a 
purely theoretical analysis. The DFs can also serve to gen- 
erate initial conditions with Monte Carlo simulators (e.g. 
iKazantzidis et al.ll2004 iBuvle et al.ll2007bl) to investigate the 
physical processes within these equilibrium models by means 
of controlled numerical A^-body simulations. Finally, our al- 
gorithm and base functions can be used in other dynami- 
cal studies, using different data moments than the density. 
This could lead to a significant improvement in the dynam- 
ical modeling of observed stellar systems, such as clusters of 
galaxies. 




APPENDIX 

DERIVATION OF THE AUGMENTED MOMENTS 

In a spherical dynamical model, the anisotropic velocity moments of the DF are 
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M2«,2m(r) = 2^M,o, / dv, / F{E,L) v"^^^ dv^ . 



Using the augmented density formalism, these moments can be calculated as [iinimif) = f^2n.2in('4'(r),r), with 
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where Df denotes the mth differentiation with respect to x. Applying this to our base functions (l44l i. we obtain for the augmented 
second-order moments 
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with y = (ip/ ipoY', and 

P'02M, r) = 2 (1 - /3(r)) jlioM.r) 
The augmented fourth-order moments have the form 
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mM^r) = 2(1- Pir)) fl4o.i(tp, r), 



Ao4,/(V','")=^ 



^ 1 - /3(r)) (2 - /3(r)) - — (/3o - /3(r)) (/3oo - (3{r)) 



The true velocity moments are connected with the anisotropic velocity moments through the relation 
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so that 



A2/,2m,2n(V','")= -B l^'" + 2 ' " 2 j /^2/,2(m+n)(V', 
N 

paj(r) = p2Qo(r) = ^ fl/v,/ /i20,i {i'(r),r) , 
i=i 

1 

paj(r) = pQ2o(r) = 2 X! ^"2,; , 

!=1 

P(^r>('") = M400('-) = ^flA',;A40,; (?/'(''),'") , 

1=1 
3 ^ 

P{vt){r) = Pimir) = -^aN.iPMj {}p{r),r) . 
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Finally, the fourth-order anisotropy profile can be derived by combining eqs. (IA7l i. ( lAl lb and ( IA12l i. 

/34W= l--|^(r)= i/3(r)(3-/3(r)) + ^ (/3o-/3(r)) {P^-l3(r)) . (A13) 
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